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Algorithms for the computation of geodesies on an ellipsoid of revolution are given. These provide accurate, 
robust, and fast solutions to the direct and inverse geodesic problems and they allow differential and integral 
properties of geodesies to be computed. 
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1. INTRODUCTION 

The shortest path between two points on the earth, custom- 
arily treated as an ellipsoid of revolution, is called a geode- 
sic. Two geodesic problems are usually considered: the direct 
problem of finding the end point of a geodesic given its start- 
ing point, initial azimuth, and length; and the inverse problem 
of finding the shortest path between two given points. Refer- 
ring to Fig.Q] it can be seen that each problem is equivalent to 
solving the geodesic triangle NAB given two sides and their 
included angle (the azimuth at the first point, a%, in the case 
of the direct problem and the longitude difference, A 12, in the 
case of the inverse problem). The framework f o r solving these 
problems w as laid down by Eegendrel (fl806b. JOrianil (Il806l 
1 18081 11810b iBessel d 18251). a nd iHelmertl dl880l) . Based on 
these works. Ivincentv ( 1975al) devised algorithms for solving 
the geodesic problems suitable for early programmable desk 




FIG. 1 The ellipsoidal triangle NAB. N is the north pole, NAF 
and NBH are meridians, and AB is a geodesic of length S12. The 
longitude of B relative to A is A12; the latitudes of A and B are (f>i 
and 4>2- EFH is the equator with E also lying on the extension of the 
geodesic AB; and arj, ati, and Q2 are the azimuths (in the forward 
direction) of the geodesic at E, A, and B. 
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calculators; these algorithms are in widespread use today. A 
good summary of Vin centy' s algor ithms and the earlier work 
in the field is given by Rapp ( 1993, Chap. 1). 

The goal of thi s paper is to adapt the geodesic methods 
of IHelmertl dl880f) and his predecessors to modern comput- 
ers. The current work goes beyond Vincenty in three ways: 
(1) The accuracy is increased to match the standard preci- 
sion of most computers. This is a relatively straightforward 
task of retaining sufficient terms in the series expansions and 
can be achieved at little computational cost. (2) A solution 
of the inverse problem is given which converges for all pairs 
of points. (Vincenty's method fails to converge for nearly an- 
tipodal points.) (3) Differential and integral properties of the 
geodesies are computed. The differential properties allow the 
behavior of nearby geodesies to be determined, which enables 
the scales of geodesic projections to be computed without re- 
sorting to numerical differentiation; crucially, one of the dif- 
ferential quantities is also used in the solution of the inverse 
problem. The integral properties provide a method for find- 
ing the area of a geodesic polygon, extending the work of 
lDanielsenldl989b . 

Section |2] reviews the classical solution of geodesic prob- 
lem by means of the auxiliary sphere and provides expan- 
sions of the resulting integrals accurate to 0(f 6 ) (where / 
is the flattening of the ellipsoid). These expansions can be 
inserted into the solution for the dire ct geodesic problem pre- 
sented by, for example, Rapp d 19931) to provide accuracy to 
machine precision. Section[3]gives the differen tial pr operties 
of geodesies reviewing the results of Helmert (1880) for the 
reduced length and geodesic scale and give the key properties 
of these quantities and appropriate series expansions to allow 
them to be calculated accurately. Knowledge of the reduced 
length enables the solution of the inverse problem by New- 
ton's method which is described in Sect. [4] Newton's method 
requires a good starting guess and, in the case of nearly antipo- 
dal points, this is pr ovided by an ap proximate solution of the 
inverse problem by Helmert (1880), as given in Sect. [5] The 
computation of area between a geodesic and t he equator is for - 
mulated in Sect. [6] extending the work of iDanielsenl (11989b . 
Some details of the implementation and present accuracy and 
timing data are discussed in Sect. [7] As an illustration of the 
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use of these algorithms, Sect. [8]gives an ellipsoidal gnomonic 
projection in which geodesies are very nearly straight. This 
provides a convenient way of solving several geodesic prob- 
lems. 

For the purposes of this paper, it is useful to generalize the 
definition of a geodesic. The geodesic curvature, n, of an arbi- 
trary curve at a point P on a surface is defined as the curvature 
of the projection of the curve onto a plane tangent to the sur- 
face at P. All shortest paths on a surface are straight, defined 
as k = at every point on the path. In the rest of this paper, 
I use straightness as the defining property of geodesies; this 
allows geodesic lines to be extended indefinitely (beyond the 
point at which they cease to be shortest paths). 

Several of the results repor ted here appeared earlier in a 
technical report, lKarnevl(l201lh . 

2. BASIC EQUATIONS AND DIRECT PROBLEM 
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FIG. 2 The elementary ellipsoidal triangle NEP mapped to the aux- 
iliary sphere. NE and NPG are meridians; EG is the equator; and 
EP is the great circle (i.e., the geodesic). The corresponding el- 
lipsoidal variables are shown in parentheses. Here P represents an 
arbitrary point on the geodesic EAB in Fig.Q] 



I consider an ellipsoid of revolution with equatorial radius 
a, and polar semi-axis b, flattening /, third flattening n, ec- 
centricity e, and second eccentricity e' given by 



/ = (a - b)/a = 1 - Vl-e 2 , 
n=(a-b)/(a + b) = f/(2-f), 
e 2 = (a 2 -b 2 )/a 2 = f(2~f), 
e' 2 = (a 2 - b 2 )/b 2 = e 2 /(l - e 2 ). 



(1) 
(2) 
(3) 
(4) 



As a consequence of the rotational symmetry of th e ellipsoid, 
geodesies obey a relation found by lClairau 3 d 17351) . namely 



sin «o = sin a\ cos /?i = sin cos /?2 , 



(5) 



where /3 is the reduced latitude (sometimes called the para- 
metric latitude), given by 



tan/3 = (! — /) tan< 



(6) 



The geodesic problems are most easily solved by using an 
auxiliary sphere which allows an exact correspondence to be 
made between a geodesic and a great circle on a sphere. On 
the sphere, the latitude <fi is replaced by the reduced latitude 
/3, and azimuths a are preserved. From Fig. it is clear that 
Clairaut's equation, sinao = sin a cos (3, is just the sine rule 
applied to the sides NE and NP of the triangle NEP and 
their opposite angles. The third side, the spherical arc length 
cr, and its opposite angle, the spherical longitude u>, are related 
to the equivalen t quantities o n the ellipsoid, the distance s and 
longitude A, by (Rappl ll993L Eqs. (1.28) and (1.170)) 



/ Vl + k 2 sin 2 a 1 da' = h(a), 
Jo 

2-/ 



(7) 



A = uj — f sin do 



■ sin 2 a' 



■.da' 



io 1 + (1 - /)VT+fc 2 sii 
ui - /sin a I 3 (cr), (8) 



where 



See also Eqs. (5.4.9) and (5.8.8) o flHelmerll(ll880l) . The origin 
for s, a, A, and uj is the point E, at which the geodesic crosses 
the equator in the northward direction, with azimuth a . The 
point P can stand for either end of the geodesic AB in Fig.Q] 
with the quantities f3, a, a, ui, s, and A acquiring a subscript 
1 or 2. I also define s\2 — S2 — si as the length of AB, with 
A12, cti2, and u>i2 defined similarly. (In this paper, «2 is the 
forward azimuth at B. Several authors use the back azimuth 
instead; this is given by c*2 ± 7r.) 

Because Eqs. ^ and ([H) depend on ctQ, the mapping be- 
tween the ellipsoid and the auxiliary sphere is not a global 
mapping of one surface to another; rather the auxiliary sphere 
should merely be regarded as a useful mathematical technique 
for solving geodesic problems. Similarly, because the ori- 
gin for A depends on the geodesic, only longitude differences, 
e.g., A12, should be used in converting between longitudes rel- 
ative to the prime meridian and A. 

In solving the spherical trigonometrical problems, the fol- 
lowing equations relating the sides and angles of NEP are 



useful, 






a = 


ph(|cosa + isinasin/3| + isinacos/3), 


(10) 


a = 


ph(cos a cos P + i sin (3) , 


(11) 


UJ = 


ph(cos cr + i sin ao sin cr) , 


(12) 


P = 


ph(|cos ao cos cr + i sinao| + i cos ao sin a), 


(13) 


a = 


ph(cos ao cos a + i sin ao ) , 


(14) 



where i = V— 1 and phfz + iy) is the phase of a complex 
number dOlver et al. i I2Q10L §1.9(1)), typically given by the 
library function atan2(y, x). Equation ([Toi l merely recasts 
Eq. (O in a form that allows it to be evaluated accurately when 
ao is close to ^71". The other relations are obtained by applying 
Napier's rules of circular parts to NEP. 

The distance integral, Eq. (0, can be expanded in a Fourier 
series 



h (cr) = Ai (a + ^2 C u sin 2 ^cr 



e cosao- 



(9) 



(15) 
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with the coefficients determined by expand ing the integra l in 
a T aylor series. It is advantageous to follow Besse| d 1 825l §5) 
and|Helmerj£l880, Eq. (5.5.1)) and use e, defined by 



Vl + fc2- i 

e = , or k 



(16) 



as the expansion parameter. This leads to expansions with 
half as many terms as the corresponding ones in fc 2 . The ex- 
pansion can be conveniently carr ied out to arbitr ary order by a 
computer algebra system such as Maxim 342009b which yields 



Ai = 


(l_ e) -l(l+l e 2 + 


64 c 


C n = 


2 ' 16 32 c ' 




C\2 = 


_J_ e 2 , J_ £ 4 _ _9_ 

16 c T 32 c 2048 


e 6 + 


Cia = 


48 ^256 ~ 




C14 = 


512 c T 512 c ~ 


J 


Cl5 = 


^e 5 + ... 

1280 ' ' 




Gl6 = 


2048 ~ 





256 c 



(17) 



(18) 



This extends Eq. (5.5.7) of|Hehriert (1880) to higher order. 
These coefficients may be inserted into Eq. (1.40) of Rapp 
dl993l) using 



Mi, for^O, 
\2Axdu for j = 21, with l> 0, 
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where here, and subsequently in Eqs. ( 1221 and ( 1261 1. a script 
letter, e.g., B, is used to stand for Rapp's coefficients. 

In the course of solving the direct geodesic problem (where 
S12 is given), it is necessary to determine a given s. Vin- 
centy solves for a iteratively. However, it is simpler to follow 
Helmertl Jl880l §5.6) and substitute s = bAiT into Eqs. 
and dTBT ), to obtain r = a + Cn sin 2la; this may be in- 
verted, for example, using Lagrange reversion, to give 



00 

E 

{=1 



C' u sin 21t, 



(20) 



where 



C' n = 


l P 9 3 1 205 5 1 
2 fc 32 fc <~ 1536 ' 


C\2 — 


5 2 37 4 , 1335 6 , 
16 96 ~ r 4096 ' 


Cl3 — 


29 3 75 5 , 

96 128 ' ' 


C 14 = 


539 4 2391 6 , 

1536 fc 2560 t ' ' ' ' 


c[ 5 = 


3467 5 , 

7680 fc ' ' 


C 16 = 


38081 6 
61440 fc ' 



(21) 



This extends Eq. (5.6.8) of lHelmertl dl880h to higher order . 
These coefficients may be used in Eq. (1.142) o fRapp(tl993|) 
using 



for j = 21, with I > 0. 



(22) 



TABLE 1 The parameters for the WGS84 ellipsoid used in the ex- 
amples. The column labeled "Eq." lists the equations used to com- 
pute the corresponding quantities. 



Qty. 


Value 


Eq. 


a 


6 378 137 m 


given 


f 


1/298.257 223 563 


given 


b 


6 356 752.314 245 m 


CD 


c 


6 371007.180 918 m 


(60} 


n 


0.001679 220 386 383 70 


£9 


e 2 


0.006 694 379 990 14132 


(3) 


e' 2 


0.006 739 496 742 276 43 





Similarly, the integral appearing in the longitude equation, 
Eq. ([§), can be written as a Fourier series 

oo 

h (a) = A 3 (a + C M sin 2laj . (23) 

(=i 



Following iHelmer 



, I expand jointly in n and e, both 



of which are O( f), to give 
i 4 3 = l_(i_i n ) e _(i + i, 



^_ e 5 

128 



(24) 



l„2^2 , ,3 , 3 n _l n 2 )e 3 



<2n = (| - \n)e + (| - \n 2 )e 2 + (£ + f; „ (> , 

~ V128 ' 64 'V e ' 128 e ' ' 



C32 

C33 

C34 
C35 



3 n + 4 n 2 le9+f#r _^ n _^ na)e 3 



JL _ _ 

,16 32'" 1 32'" 1 V64 32 

+ (lis + T28 n ) € + 2§6 £ + '"> 



V 192 64"^ 192" / c ^ V128 192 ul 

+ -J- f 5 J 

T 512 fc ~ ' 

V512 256 V c ' 512 c ~ 



21 5 
2560 



(25) 



This extends Eq. (5.8.14) of lHelmertl dl880h to higher order 
These coefficients may be inserted into Eq. (1.56) of iRappI 
d 19931) using 



Aj — 




forj = 0, 

for j = 21, with / > 0. 



(26) 



The equations given in this section allow the direct geode- 
sic problem to be solved. Given fa (and hence fa) and ax 
solve the spherical triangle NEA to give ao, o\, and oj\ using 
Eqs. ([Tol l. (Hit , and (fT2] i. Find si and Ai from Eqs. (O and 
© together with Eqs. ( fT3T > and ( |23l . (Recall that the origin 
for A is E 1 in Fig. [T]) Determine sq, = Si + S12 and hence 
(j 2 using Eq. d20b . Now solve the spherical triangle NEB to 
give ct2, fa (and hence fa), and k>2, using Eqs. (fl~4l . dT3l . and 
(fl~2] i. Finally, determine A2 (and A12) from Eqs. ([8]l and d23l . 
A numerical example of the solution of the direct problem is 
given in Table fusing the parameters of Table Q] 
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TABLE 2 A sample direct calculation specified by 4>\ = 40°, Qi = 
30° , and S12 = 10 000 km. For equatorial geodesies (<j)i = and 
ot\ = 2 7r )> Eq< tH) i s indeterminate; in this case, take o\ = 0. 



Qty. 


Value 


Eq. 


<t>i 


40° 


given 


Ql 


30° 


given 


S12 


10 000 000 m 


given 


Solve triangle NEA 




Pi 


39.905 277146 01° 


ID 


Oto 


22.553 940 202 62° 




<Tl 


43.999 153 645 00° 


jiii 

qjj 


Ul 


20.323 718 278 37° 


CCD 


Determine <T2 




k 2 


0.005 748 029 628 57 


(2) 


c 


0.001 432 892 204 16 


EB 


Ai 


1.001435 462 362 07 


CCD 


h{ai) 


0.768 315 388 86412 




Si 


4883 990.626 232 m 


CD 


S2 


14 883 990.626 232 m 


Sl + S12 


T2 


133.962 660 502 08° 


s 2 /{bAi) 


C"2 


133.921 640 830 38° 


CIS 


Solve triangle NEB 




a 2 


149.090 169 318 07° 


1 , , 1 

11141 


/3a 


41.697 718 092 50° 


CCD 


W2 


158.284 121471 12° 


CCD 


Determine A12 




A3 


0.999 284 243 06 


en 


/ 3 (<Tl) 


0.767 737 860 69 


CSJ 


J r 3(o"2) 


2.335 343 22170 


ED 


Ai 


20.267 150 380 16° 


cd 


A 2 


158.112 050 423 93° 


ID 


A12 


137.844 900 043 77° 


A2 — Ai 


Solution 




02 


41.793 310 205 06° 


(6) 


A12 


137.844 900 043 77° 




a 2 


149.090 169 318 07° 





3. DIFFERENTIAL QUANTITIES 

Before turning to the inverse problem, I present Gauss' so- 
lution for the differential behavior of geodesies. One differen- 
tial quantity, the reduced length m\i, is needed in the solution 
of the inverse problem by Newton's method (Sect. HI and an 
expression for this quantity is given at the end of this section. 
However, because this and other differential quantities aid in 
the solution of many geodesic problems, I also discuss their 
derivation and present some of their properties. 

Consider a reference geodesic parametrized by distance s 
and a nearby geodesic s eparated from the reference by in- 
finitesimal distance t(s). iGaussI (Q828) showed that t(s) sat- 
isfies the differential equation 



dH(s) 
ds 2 




FIG. 3 The definitions of mvz and M12 are illustrated in (a) and (b). 
A geometric proof of Eq. J29l > is shown in (c); here AB and A'B' 
are parallel at B and B', BAB' = dai, BB' = m 12 dai, AA' = 
M2imi2 dai, and finally AB' A' — M21 dai, from which Eq. ([29} 
follows. 



where K(s) is the Gaussian curvature of the surface. As a 
second order, linear, homogeneous differential equation, its 
solution can be written as 

t(s) = At A (s) + Bt B (s), 

where A and B are (infinitesimal) constants and Ia and ts are 
independent solutions. When considering the geodesic seg- 
ment spanning si < s < S2, it is convenient to specify 



t A (si) = 0, 
*b(si) = 1, 



dt A (s) 



ds 
dt B (s) 



d.s 



= 1, 



0, 



and to write 



The quantity mi 9 is the reduced length of the geodesic 
(IChristoffelL [[868). Consider two geodesies which cross at 
s = si at a small angle dai, Fig[3ja); at s = s 2 , they will 
be separated by a distance mi 2 dai. Similarly I call Mi 2 the 
geodesic scale. Consider two geodesies which are parallel at 
s = si and separated by a small distance dti, Fig[3jb); at 
s — s 2 , they will be separated by a distance M\ 2 dt\. 

Several relations between nii 2 and Mi 2 follow from the 
defining equation, Eq. ( |27|). The re duced length obeys a reci- 
procity relation JChristoffell[l868l §9), m 2 i + TO12 = 0; the 
Wronskian is given by 



W(M 12 ,mi 2 )(s 2 ) = Mi 



dmi 2 

> 

ds 2 



dM 12 
TO12— : — = 1; 
ds 2 



+ K(s)t(s)=0, 



(27) 



and the derivatives are 

dm i2 

ds 2 
dM 12 
ds 2 



= M 2U 

1 - Mi 2 M 2 i 
m 12 



(28) 

(29) 
(30) 
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The constancy of the Wronskian follows by noting that its 
derivative with respect to s 2 vanishes; its value is found by 
evaluating it at s 2 — Si. A geometric proof of Eq. (l29l i is 
given in Fig |3jc) and Eq. (f30b then follows from Eq. d28l) . 
With knowledge of the derivatives, addition rules for m,\i and 
M12 are easily found, 



M12 = COS tTi COS (72 



\/\ + k 2 si" 2 



\/l + fc 2 sin 2 (Ti 



sin (72 

sin cr i sin (72 



Sin(7l COS(7 2 (j((72) - J((7l)) 

\pY + k 2 sin 2 ai 



(39) 



where 



7™13 


= m 12 M 2 3 


+ BI23M21, 


(31) 




M 13 


= M 12 M 23 


-(1-M 12 M 21 ) — , 


(32) 


J((7) = f 
Jo 


M 31 


= M 32 M 21 


-(I-M23M32) — , 
m 2 3 


(33) 


s 

= b~ 



k 2 sin 2 a' 



.da' 



■.da' 



where points 1, 2, and 3 all lie on the same geodesic. 

Geodesies allow concepts from plane geometry to be gen- 
eralized to apply to a curved surface. In particular, a geodesic 
circle may be defined as the curve which is a constant geode- 
sic distance from a fixed point. Similarly, a geodesic parallel 
to a reference curve is the curve which is a constant geodesic 
distance from that curve. (Thus a circle is a special case of a 
parallel obtained in the limit when the reference curve degen- 
erates to a point.) Parallels occur naturally when considering, 
for example, the "12-mile limit" for territorial waters which 
is the boundary of points lying within 12 nautical miles of a 
coastal state. 

The geodesic curvature of a parallel can be expressed in 
terms of m\ 2 and M\ 2 . Let point 1 be an arbitrary point on 
the reference curve with geodesic curvature k\. Point 2 is the 
corresponding point on the parallel, a fixed distance s\ 2 away. 
The geodesic curvature of the parallel at that point is found 
from Eqs. (|29]i and (O, 



K 2 



M 2lKl - (1 - M 12 M 21 )/m 12 
mi 2 Ki + M12 



(34) 



The curvature of a circle is given by the limit k,\ — > 00, 

K 2 = M 21 /m 12 . (35) 

If the reference curve is a geodesic (ki —> 0), then the curva- 
ture of its parallel is 



k 2 = -(1 - M 12 M 21 )/{M 12 m 12 ). 



(36) 



If the reference curve is indented, then the parallel intersects 
itself at a sufficiently large distance from the reference curve. 
This begins to happen when n 2 — > 00 in Eq. (|34l . 

The results above apply to general surfaces. For a geodesic 
on an ellipsoid of revolution, the Gaussian curvature of the 
surface is given by 



K = 



(1 — e sin 



2 r.;„2 i\2 



1 



(37) 



b 2 6 2 (l + fc 2 sin 2 (7) 2 ' 

iHelmertld 18801 Eq. (6.5.1)) solves Eq. (l27l i in this case to give 



i\ 2 /b = V 1 + k 2 sin 2 (7 2 coscti sin 02 



\/ \ + k 2 sin 2 (7i sin o\ cos a 2 

C0S(7i C0S(72(>/(C2) — ^(ci))j 



(38) 



h(a)-I 2 {a). 



(40) 



Equati on ( |39l may be obtained from Eq. (6.9.7) of iHelmertl 
(1880), which gives dmi 2 /ds 2 ; M\ 2 may then be found from 
Eq. ( 1291 with an interchange of indices. In the spherical limit, 
/ ->■ 0, Eqs. (J38]> and ([39]) reduce to 

TO12 = asin<7i2 = asin(si2/a), 
M i2 = cos CT12 = cos(si 2 /a). 

The integral I 2 (a) in Eq. d40l may be expanded in a Fourier 
series in similar fashion to ii(<7), Eq. (1151 l. 

Ha) =A 2 (<J + J2 ° 21 sin 2?fT )' (41) 



z=i 



where 



A 2 = (l- £ )(l + I e 2 + JL e 4 + iL e 6 + ...) I (42) 

^21 - 2 £ + IE 6 + 32 £ H ' 

C = 3-c 2 A- -Le 4 -J- 35 r 6 4. 

^22 16 c T 32 c T 2048 ~ ' 

L-23 — TS- £ + ^H7j£ 

C*24 = 

^ 2S = lli.su 
C26 = 



48" 1 256" 

512 ' 512 c 
63 £ 5 + . 



77 6 
2048 - 



(43) 



4. INVERSE PROBLEM 



The inverse problem is intrinsically more complicated than 
the direct problem because the given included angle, A12 in 
Fig. Q] is related to the corresponding angle on the auxiliary 
sphere oj\ 2 via an unknown equatorial azimuth cto. Thus, the 
inverse problem inevitably becomes a root-finding exercise. 

I tackle this problem as follows. Assume that a\ is known. 
Solve the hybrid geodesic problem: given <fii, <fi 2 , and ax, find 
A12 corresponding to the first intersection of the geodesic with 
the circle of latitude (j> 2 . The resulting A12 differs, in general, 
from the given Ai 2 ; so adjust a\ using Newton's method until 
the correct A12 is obtained. 

I begin by putting the points in a canonical configuration, 



h <o, 



< 



< A12 < 77. (44) 
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a, (°) 



FIG. 4 The variation of A12 as a function of ai for (f>\ = —30°, 
various fc, and the WGS84 ellipsoid. Part (a) shows A12 for <j>2 = 
0°, ±15°, ±25°, and ±30°. For \cj> 2 \ < \<j>i\, the curves are strictly 
increasing, while for 4>2 = ±01, the curves are non-decreasing with 
discontinuities in the slopes at ai = 90°. An enlargement of the top 
right corner of (a) is shown in (b) with <j>2 G [29°, 30°] at intervals 
of 0.2°. 



This may be accomplished swapping the end points and the 
signs of the coordinates if necessary, and the solution may 
similarly be transformed to apply to the original points. All 
geodesies with a,\ € [0,7r] intersect latitude <fr 2 with A12 G 
[0, 7r]. Furthermore, the search for solutions can be restricted 
to a 2 6 [0, |7r], because this corresponds to the first intersec- 
tion with latitude (j> 2 ■ 

Meridional (A12 = or tt) and equatorial (<f>i = cf> 2 = 0, 
with A12 < (1 — f)if) geodesies are treated as special cases, 
since the azimuth is then known: ai — A12 and ai = 
respectively. The general case is solved by Newton's method 
as outlined above. 

The solution of the hybrid geodesic problem is straightfor- 
ward. Find (3i and (3 2 from Eq. ©, solve for ao and a 2 from 
Eq. (0, taking cos olq > and cos a 2 > 0. In order to com- 
pute a 2 accurately, use 



cosa2 



+ \/cos 2 a\ cos 2 j3i + (cos 2 j3 2 — cos 2 Pi 
cos 



(45) 



in addition to Eq. ((5]). Compute ai, 0J1, and 0J2 using 
Eqs. ( fTTT ) and ( fT2b . Finally, determine A12 (and, once conver- 
gence is achieved, S12) as in the solution to the direct problem. 
The behavior of A12 as a function of ai is shown in Fig. [4] 

To apply Newton's method, an expression for dAi2/dai is 
needed. Consider a geodesic with initial azimuth ai. If the 
azimuth is increased to ai + dai with its length held fixed, 
then the other end of the geodesic moves by 77112 dai in a 
direction ^tt + a 2 - If the geodesic is extended to intersect 
the parallel fa once more, the point of intersection moves by 
m^dai/ cosa2; see Fig- The radius of this parallel is 
a cos $2 ; thus the rate of change of the longitude difference is 



dAi 



mi2 



1 



dai 



a cos a 2 cos (3 2 



(46) 



This equation can also be obtained from Eq. (6.9.8b) of 
iHelmertl ([T880). Equation d46*l > becomes indeterminate when 
/3 2 = ±/3i and ai — because mia and cosa2 both van- 
ish. In this case, it is necessary to let ai 
the limit S —> ±0, which gives 



t7t + 5 and to take 



dAi 



dai 



\J\ — e 2 cos 2 pi 
sin (3i 



[1 Tsign(cosai)), (47) 



where sign(cosai) = — sign(<5). A numerical example of 
solving the inverse geodesic problem by this method is given 
at the end of the nex t section. 

I Vincentvl d 1 975ah . who uses the iterative method of lHelmerj 
(118801 S5 .13) to solve the inverse problem, was aware of its 
failure to con verge for nearly a ntipodal points. In an unpub- 
lished report dVincentvi Il975bl) . he gives a modification of 
his method which deals with this case. Unfortunately, this 
sometimes requires many thousands of iterations to converge, 
whereas Newton's method as described here only requires a 
few iterations. 



5. STARTING POINT FOR NEWTON'S METHOD 

To complete the solution of the inverse problem a good 
starting guess for ai is needed. In most cases, this is provided 
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TABLE 3 First sample inverse calculation specified by cj>\ = 
-30.123 45°, fa = -30.123 44°, and \ 12 = 0.000 05°. Because 
the points are not nearly antipodal, an initial guess for a\ is found as- 
suming oji2 = A12 /w. However, in this case, the line is short enough 
that the error in CJ12 is negligible at the precision given and the so- 
lution of the inverse problem is completed by using S12 = awai2. 
More generally, the value of ai would be refined using Newton's 
method. 



Qty. 


Value 


h,q. 


fa 


-30.123 45° 


given 


fa 


-30.123 44° 


given 




0.000 05° 


given 


Determine W12 




Pi 


-30.039 990 838 21° 


(6) 


02 


-30.039 980 854 91° 


© 


W 


0.997 488 477 44 


idiii 


Wl2 


0.000 050125 89° 


Aia/io 


<T12 


0.000 044 526 41° 


ED 


Solution 






Ol 


77.043 533 542 37° 


@9li 




77.043 508 449 13° 


EB 


S12 


4.944 208 m 





by assuming that CJ12 = X12/1V, where 



w = \Jl- e 2 ((cos/3i +cos/? 2 )/2)' 



(48) 



and solving for th e great circle on the auxiliary sphere, using 
dVincentylll975ah 



Z\ = cos fa sin fa — sin fa cos fa cos W12 

+ icos fa sina;i2, 
Z2 = — sin fa cos fa + cos fa sin fa cos W12 

+ i cos fa sina;i2, 
at = phzi, 
a- 2 = phz 2 , 

<j\2 = ph(sin fa sin fa + cos fa cos /?2 cos U)t2 + i |#i 



(49) 
(50) 
)■ 

(51) 



An example of the solution of the inverse problem by this 
method is given in Table [3] 

This procedure is inadequate for nearly antipodal points be- 
cause both the real and imaginary components of z\ are small 
and at depends very sensitively on W12. In the correspond- 
ing situation on the sphere, it is possible to determine ot\ by 
noting that all great circles emanating from A meet at O, the 
point antipodal to A. Thus at may be determined as the sup- 
plement of the azimuth of the great circle BO at O; in addi- 
tion, because B and O are close, it is possible to approximate 
the sphere, locally, as a plane. 

The situation for an ellipsoid is slightly different because 
the geodesies emanating from A, instead of meeting at a point, 
form an envelope, centered at O, in the shape of an astroid 





a 
















-x/(1+u) 



-yl\i 



FIG. 6 The solution of the astroid equations by similar triangles. 
The scaled coordinates of B are (x,y); O is the point antipodal to 
A. The line BCD, which is given by Eq. d54t . is the continuation of 
the geodesic from AB with C being its intersection with the circle 
f3 — —fa and D its intersection with the meridian A = Ai + n. The 
envelope of lines satisfying CD — 1 gives the astroid, a portion of 
which is shown by the curves. 



whose extent is O(f) d Jacobill 1 89 it Eqs. (16)-(17)). The po- 
sition at which a particular geodesic touches this envelope is 
given by the condition mi2 = 0. However elementary meth- 
ods can be used to determine the envelope. Consider a geode- 
sic leaving A (with fa < 0) with azimuth ct\ G [hn, n]. This 
first intersects the circle of opposite latitude, fa = — fa, with 
er 12 = W12 = 7T and a 2 = tt — a.\. Equation ([H} then gives 



A12 = 7r — fn cos fa sin a\ + 0(/ 2 



(52) 



Define a plane coordinate system (x, y) centered on the an- 
tipodal point where A = fair cos 2 fa is the unit of length, 
i.e., 



A12 



a cos fa 



fa 



-fa 



-y- 



(53) 



In this coordinate system, Eq. d52l i corresponds to the point 
x = — sin at, y = and the slope of the geodesic is — cot at. 
Thus, in the neighborhood of the antipodal point, the geodesic 
may be approximated by 



y 



smai cosai 



1 = 0, 



(54) 



where terms of order / 2 have been neglected. Allowing at 
to vary, Eq. d54l ) defines a family of lines approximating the 
geodesies emanating from A. Differentiating this equation 
with respect to at and solving the resulting pair of equations 
for x and y gives the parametric equations for the astroid, 
x = — sin 3 ai, and y — — cos 3 ai. Note that, for the or- 
dering of points given by Eq. d44K x < and y < 0. 

Given x and y (i.e., the position of point B), Eq. (|54| ) may 
be solved to obt ain a first appr oximation to at- This prescrip- 
tion is given bv lHelmertl (fl 880, Eq. (7.3.7)) who notes that this 
results in a quartic which may be found using the construction 
given in Fig. [6] Here COD and BED are similar triangles; 
if the (signed) length BC is \i, then an equation for /1 can be 
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TABLE 4 Second sample inverse calculation specified by 0i = 
-30°, 02 = 29.9°, and Ai 2 = 179.8°. Because the points are 
nearly antipodal, an initial guess for qi is found by solving the as- 
troid problem. Here /i is the positive root of Eq. i55\ , If y — 0, then 
ai is given by Eq. i57i . The value of ct\ is used in Table[3] 



TABLE 5 Second sample inverse calculation, continued. Here 



Qty. 


Value 


Eq. 


01 


-30° 


given 


02 


29.9° 


given 


Al2 


179.8° 


given 


Solve the astroid problem 


X 


-0.382 344 




V 


-0.220 189 




A 4 


0.231 633 


m 


Initial 


guess for qi 




Ql 


161.914° 





found by applying Pythagoras' theorem to COD, 



= 1, 



which can be expanded to give a 4th-order polynomial in /.i, 

+ 2 M 3 + (1 - x 2 - y V - 2y 2 n -y 2 =0. (55) 

Descartes' ru le of signs shows that, for y ^ 0, there is one 
positive root (lOlver et a/.ll20ld §l.ll(ii )) and this is the solu- 
tion corresponding to t he shortest path. T his root can be found 
by standard methods (lOlver et ali 1201 Oi §l.ll(ii i)). Equa- 
tion d55l l arises in converting from geocentric to geodetic co- 
ordinates, and I use the solution to that problem given by 
IVermeillg d2002l) . The azimuth can then be determined from 
the triangle COD in Fig.|6] 

a x = ph(y/n - ix/(l + ^)). (56) 

If y = 0, the solution is found by taking the limit y — >• 0, 

«i = ph(±-\/ max(0, 1 — x 2 ) — ix) . (57) 

Tables [4^6] together illustrate the complete solution of the in- 
verse problem for nearly antipodal points. 



6. AREA 

The last geodesic algorith m I presen t is for geodesic areas. 
Here, I extend the method of Danielsen ( 1989) to higher order 
so that the result is accurate to round-off, and I recast his series 
into a simple trigonometric sum. 

Let S\2 be the area of the geodes ic quadrilateral AFHB in 
Fig.Q] Following|Darrielsen|(l989), this can be expressed as 
the sum of a spherical term and an integral giving the ellip- 
soidal correction, 



\\°2 denotes the desired value of the longitude difference; Newton's 



method is used to adjust a\ so that A12 
ai is used in Table|6] 



A 



(0) 



The final value of 



Qty. 


Value 


Eq. 


0i 


-30° 


given 


02 


29.9° 


given 


ai 


161.914° 


Tableg] 


\(0) 
A 12 


179.8° 


given 


Solve triangle NEA 




fa 


-29.916 747 713 24° 


© 


Cl() 


15.609 39746414° 


GoJ 




-148.812 535 665 96° 


{TQ 


Ul 


-170.748 966 96128° 




Solve trianj 


;le NEB 




ft 


29.816 916 42189° 


© 


a 2 


18.067 287962 31° 


EKES 


(72 


31.082 449 768 95° 


(D) 


UJ2 


9.213 457611 10° 


C3 


Determine A12 




A: 2 


0.006 251537 916 62 




c 


0.001558 018 267 80 




Ai 


-170.614 835 524 58° 


m 


A 2 


9.185 420 098 39° 




A12 


179.800 255 622 97° 


A2 — Ai 


Update Qi 




\ \ (°) 

A12 — A i2 


8\i2 


0.000 255 622 97° 




-0.009 480 409 276 40 


gOji 


J(a 2 ) 


0.000 313 491 286 30 


gob 


mi2 


57 288.000 110 m 




dAi2/dai 


0.010 889 317161 15 


(O 


Sai 


-0.023 474 655 19° 


-5Ai 2 /(dAi 2 /dQi) 


Ql 


161.890 525 344 81° 


Qi + Sai 


Next iteration 




5\l2 


0.000 000 006 63° 




Ctl 


161.890 524 736 33° 





S(a) — c a + e a cosaosinao/4(cr), (59) 



where 



a 2 b 2 tanh 1 



2 2 e 



(60) 



is the authalic radius, 



r t{e' 2 )-t{k 2 sin 2 a') sm*> 
h{a) = - I 7, ,0.9, — dcr . (61) 



tt/2 e 



12 - k 2 sin 2 a> 2 



t(x) = x + v x 1 + 1 sinh X \fx. 

Expanding the integrand in powers of e' 2 and k 2 and perform- 
ing the integral gives 



S12 = S(<j 2 ) - S(a\), 



(58) 



= J^Cu cos((2? + 1) 



(62) 
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TABLE 6 Second sample inverse calculation, concluded. Here the 
hybrid problem (<f>i, <f>2, and ai given) is solved. The computed 
value of A12 matches that given in the specification of the inverse 
problem in Table|4] 



Qty- 


Value 


Eq. 


<h 


-30° 


given 




29.9° 


given 


Oil 


161.890 524 736 33° 


Table Rl 


Solve triangle NEA 




Pi 


-29.916 747 713 24° 


,,,, 


ao 


15.629 479 665 37° 




01 


-148.809 136 917 76° 


ifTTIi 
n 1 1 V 


Wl 


-170.736 343 780 66° 


ILZl 


Solve trian 


gle NEB 




ft 


29.816 916 42189° 




Q.2 


18.090 737 245 74° 




ai 


31.085 834 470 40° 


CD 




9.226 028 621 10° 




Determine sja and A12 




Sl 


-16 539 979.064 227 m 


(7) 


S2 


3 449 853.763 383 m 





S12 


19 989 832.827 610 m 


S2 ~ Sl 


Ai 


-170.602 04712148° 


m 


A 2 


9.197 952 878 52° 




A12 


179.800 000 000 00° 


A2 — Ai 



Solution 




ai 


161.890 524 736 33° 


Oil 


18.090 737 245 74° 


S12 


19 989 832.827 610 m 



TABLE 7 The calculation of the area between the equator and the 
geodesic specified by <f>i — 40°, ct\ — 30°, and S12 = 10 000 km. 
This uses intermediate values computed in Table|2] 



Otv 

v L y- 


Value 


Fn 
nq. 


ao 


22.553 940 202 62° 


Tabled 


ai 


30° 


Tabled 


012 


149.090 169 318 07° 


Tabled 


1 


4^ qqq 1 p,a^ nn° 


TahlplTl 
1 ctuic iz.i 


(72 


133.921640 830 38° 


Table |2] 


k 2 


0.005 748 029 628 57 


Tabled 


Computi 


; area 




h{<?i) 


0.479 018 145 20 






-0.461917119 02 


462) 


S(ci) 


21298 942.66715 km 2 


m 


S(a 2 ) 


105 574 566.089 50 km 2 


m 


s 12 


84 275 623.422 35 km 2 


(EE) 



An example of the computation of S12 is given in Table|7] 

Summing S12, Eq. d58l ), over the edges of a geodesic poly- 
gon gives the area of the polygon provided that it does not 
encircle a pole; if it does, 2irc 2 should be added to the result. 
The first term in Eq. d59l contributes c 2 (a2 — cei) to Si2- This 
is the area of the quadrilateral AFHB on a sphere of radius c 
and it is proportional to its spherical excess, d2 — ocx, the sum 
of its interior angles less 2ir. It is important that this term be 
computed accurately when the edge is short (and ot\ and 012 
are nearly e qual). A suitable identity for 012 — ot\ is given by 
lBesselldl825L §11), 



where 

C40 = 



15 



4 »/4 
105 c 



315 c 



64 IS _ 128 a0\ 
3465 l " 9009 e ) 



20 



35 c 



_2_ p '4 W_J6 

105 c 1155 c 



_32_/8\ 

3003 c ) 



+ (i - J-e' 2 + -§-e' 4 - -^-e' 6 )k 4 

~ V42 63 T 693 9009 c / 



J-e' 2 

99° 



~ V110 143 



143" 

^41 = (tIo - 315 e ' 2 



_i_ 10 „/4 
~ t " 1287 

J4 16 



945 



10395' 



27027 



+ 



^252 
*360 



^-e' 2 

378 c 

^-e' 2 

495 c 



2079° 
1287 J 



4° -p /6 U 4 



27027 1 



(J, 2_ p /2\jL8 , 5 1 1 

V495 1287 c / ' 3276 



G 



12 



v495 

^2100 ~~ 

- f-L 

V1800 



~ V1925 
^43 = ( T7640 ~ 

- f^— 

V 10780 

C44 = ( ' 

a 



1287 

3150 c 

2475 c 
2 

5005 1 

24255 f 
1 



„'4 



1 17325 
12 1 2 
~T 6435 c 

jus L_£ 

2184"' 

1 2 



45045 
/4 U 6 



3 )fc 4 



,/2 



10 



12 



' 63063 
„'2\ 7,8 



e' 4 U 6 



14014 
1 



,'2\l,8 



45 



124740 162162" 
1 ;„10 , 



k< 



,10 



45864 

1 jUO 
58968 



792792 



(63) 



tan 



a 2 — o;i 



sin i(/? 2 + /3i) W12 



tan— -. (64) 



7. IMPLEMENTATION 



The algorithms described in the preceding sections can be 
readily converted into working code. The polynomial expan- 
sions, Eqs. Q3. CDU, <ED, (ED, (E3, ©, (S3, and d63]), 
are such that the final results are accurate to 0(/ 6 ) which 
means that, even for f — tI^, the truncation error is smaller 

J 150 

than the round-off error when using IEEE double precision 
arithmetic (with the fraction of the floating point number rep- 
resented by 53 bits). For speed and to minimize round-off 
errors, the polynomials should be evaluated with the Horner 
method. The parenthetical expressions in Eqs. (EH>, d25l ). and 
d63l depend only on the flattening of the ellipsoid and can be 
computed once this is known. When determining many points 
along a s ingle geodesic, t he polynomials need be evaluated 
just once. IClenshawld 19551) summation should be used to sum 
the Fourier series, Eqs. ( fT3T >, d23l , d4!l . and d62| >. 

There are several other details to be dealt with in imple- 
menting the algorithms: where to apply the two rules for 
choosing starting points for Newton's method, a slight im- 
provement to the starting guess Eq. d56l l, the convergence cri- 
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terion for Newton's method, how to minimize round-off errors 
in solving the trigonometry problems on the auxiliary sphere, 
rapidly computing intermediate points on a geodesic by using 
<7i2 as the metric, etc. I refer the reader t o the impleme nta- 
tions of the algorithms in Geographic Lib dKarnevl 1201 2|) for 
possible ways to address these issues. The C++ implemen- 
tation has been tested against a large set of geodesies for the 
WGS84 ellipsoid; this was generated by continuing the se- 
ries expansions to 0(/ 30 ) and by solving the direct problem 
using with high-precision arithmetic. The round-off errors 
in the direct and inverse methods are less than 15 nanome- 
ters and the error in the computation of the area S\2 is about 
0.1m 2 . Typically, 2 to 4 iterations of Newton's method are 
required for convergence, although in a tiny fraction of cases 
up to 16 iterations are required. No convergence failures are 
observed. With the C++ implementation compiled with the 
g++ compiler, version 4.4.4, and running on a 2.66 GHz Intel 
processor, solving the direct geodesic problem takes 0.88 /is, 
while the inverse problem takes 2.34 fis (on average). Several 
points along a geodesic can be computed at the rate of 0.37 /is 
per point. These times are comparable to those for Vincenty's 
algorithms implemented in C++ and run on the same architec- 
ture: 1.11 /us for the direct problem and 1.34 /is for the inverse 
problem. (But note that Vincenty's algorithms are less accu- 
rate than those given here and that his method for the inverse 
problem sometimes fails to converge.) 

8. ELLIPSOIDAL GNOMONIC PROJECTION 

As an application of the differential properties of geode- 
sies, I derive a generalization of the gnomonic projection to 
the ellipsoid. The gnomonic projection of the sphere has the 
property that a ll geodesies on the sphere map to straight lines 
dSnvdeii Il987l §22). Such a projection is impossible for an 
ellipso id because it d oes not have constant Gaussian curva- 
ture dBeltramil 1 1 865l §18); nevertheless, a projection can be 
constructed in which geodesies are very nearly straight. 

The spherical gnomonic projection is the limit of the dou- 
bly azimuthal projection of the sphere, wherein the bearings 
from two fixed points A and A' to B a r e pres erved, as A' 
approaches A (Bugavev skiv and Snvderl 1 1995b . The con- 
struction of the generalized gnomonic projection proceeds 
in the same way; see Fig. [7] Draw a geodesic A'B' such 
that it is parallel to the geodesic AB at A. Its initial sep- 
aration from AB is sin7df; at B 1 , the point closest to B, 
the separation becomes Mi2sin7cli (in the limit dt —> 0). 
Thus the difference in the azimuths of the geodesies A'B and 
A'B' at A' is (Af^/ro^) sin7di, which gives 7 + 7' = 
ir — (A/12/77112) sin7df. Now, solving the planar triangle 
problem with 7 and 7' as the two base angles gives the dis- 
tance AB on the projection plane as 77112/A/12. 

This leads to the following specification for the generalized 
gnomonic projection. Let the center point be A; for an arbi- 
trary point B, solve the inverse geodesic problem between A 
and B; then B projects to the point 

x = psinai, y = pcosa\, p = raYijM\2\ (65) 



B 




B' 



A 6t A' 



FIG. 7 The construction of the generalized gnomonic projection as 
the limit of a doubly azimuthal projection. 



the projection is undefined if M12 < 0. In the spherical limit, 
this becomes th e standard gnomonic projection, p = atancri2 
JSnvder[[l987l p. 165). The azimuthal scale is I/M12 and the 
radial scale, found by taking the derivative dp/dsi2 and using 
Eq. d28l ), is 1/M 2 2 . The reverse projection is found by com- 
puting ai = ph(y + ix), finding s 12 using Newton's method 
with dp/dsi2 = 1/A-f 2 2 (i-e., the radial scale), and solving the 
resulting direct geodesic problem. 

In order to gauge the usefulness of the ellipsoidal gnomonic 
projection, consider two points on the earth B and C, map 
these points to the projection, and connect them with a straight 
line in this projection. If this line is mapped back onto the 
surface of the earth, it will deviate slightly from the geodesic 
BC. To lowest order, the maximum deviation h occurs at the 
midpoint of the line segment BC; empirically, I find 

h=^(V*f-t)t, (66) 

where / is the length of the geodesic, K is the Gaussian cur- 
vature, VK is evaluated at the center of the projection A, and 
t is the perpendicular vector from the center of projection to 
the geodesic. The deviation in the azimuths at the end points 
is about Ah /I and the length is greater than the geodesic dis- 
tance by about |/i 2 /7. In the case of an ellipsoid of revolution, 
the curvature is given by differentiating Eq. (T37T > with respect 
to <fi and dividing by the meridional radius of curvature to give 

VA" = -^e 2 (l - e 2 sin 2 (/>) 5/2 cos sin </>0, (67) 

where is a unit vector pointing north. Bounding the mag- 
nitude of h, Eq. (1661 . over all the geodesies whose end points 
lie within a distance r of the center of projection, gives (in the 
limit that / and r are small) 

h f r 3 

~ < is- (68) 

r 8 a 6 

The maximum value is attained when the center of projection 
is at <fi — ±45° and the geodesic is running in an east-west 
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center of projection; iterate this process until OW — O^ -1 ' 
which is then the desired intersection point. 

The second problem is the interception problem: given a 
geodesic between A and B, find the point O on the geodesic 
which is closest to a given point C. The solution is similar 
to that for the intersection problem; however the interception 
point in the projection is 

c • (b — a)(b — a) — (a x b • z)z x (b — a) 



Provided the given points lie within about a quarter meridian 
of the intersection or interception points (so that the gnomonic 
projection is defined), these algorithms converge quadratically 
to the exact result. 



FIG. 8 The coast line of Europe and North Africa in the ellipsoidal 
gnomonic projection with center at (45°N, 12°E) near Venice. The 
graticule lines are shown at multiples of 10°. The two circles are 
centered on the projection center with (geodesic) radii of 1000 km 
and 2000 km. The data for the coast lines is taken from GMT 
dWessel and Smithll201Ch at "low" resolution. 



direction with the end points at bearings ±45° or ±135° from 
the center. 

Others have pr oposed different gene ralizations of the gno- 
monic projection. iBowringl (1 1997b and IWilliamsl ( 1 19971) give 
a projection i n whi ch great ellipses project to straight lines; 
iLetovaFtsevI ( 1 19631) suggests a projection in which normal 
sections through the center point map to straight lines. Em- 
pirically, I find that h/r is proportional to r/a and r 2 /a 2 for 
these projections. Thus, neither does as well as the projec- 
tion derived above (for which h/r is proportional to r 3 /a 3 ) at 
preserving the straightness of geodesies. 

As an illustration of the properties of the ellipsoidal gno- 
monic projection, Eq. d65l l. consider Fig.[8]in which a projec- 
tion of Europe is shown. The two circles are geodesic circles 
of radii 1000 km and 2000 km. If the geodesic between any 
two points within one of these circles is estimated by using a 
straight line on this figure, the deviation from the true geode- 
sic is less than 1.7 m and 28 m, respectively. The maximum 
errors in the end azimuths are 1.1" and 8.6" and the maximum 
errors in the lengths are only 5.4 /im and 730 fim. 

The gnomonic projection can be used to solve two geodesic 
problems accurately and rapidly. The first is the intersection 
problem: given two geodesies between A and B and between 
C and D, determine the point of intersection, O. This can be 
solved as follows. Guess an intersection point and use 
this as the center of the gnomonic projection; define a, b, c, 
d as the positions of A, B, C, D in the projection; find the 
intersection of AB and CD in the projection, i.e., 



(c x d • z)(b — a) — (a x b • z)(d — c) 
(b — a) x (d — c) • z 



(69) 



where * indicates a unit vector (a = a/a) and i — x x y is 
in the direction perpendicular to the projection plane. Project 
o back to geographic coordinates and use this as a new 



9. CONCLUSIONS 

The classical geodesic problems entail solving the ellip- 
soidal triangle NAB in Fig. Q] whose sides and angles are 
represented by fa, fa, s%2 and ct\, a-x, \vz- In the direct prob- 
lem fa, ot\, and s\2 are given, while in the inverse problem 
fa, A12, and fa are specified; and the goal in each case is to 
solve for the remaining side and angles. The algorithms given 
here provide accurate, robust, and fast solutions to these prob- 
lems; they also allow the differential and integral quantities 
mi2, M12, M21, and S12 to be computed. 

Much of the work described here involves applying stan- 
dard computational techniques to earlier work. However, at 
least two aspects are novel: (1) This paper presents the first 
complete solution to the inverse geodesic problem. (2) The 
ellipsoidal gnomonic projection is a new tool to solve various 
geometrical problems on the ellipsoid. 

Furthermore, the packaging of these various geodesic capa- 
bilities into a single library is also new. This offers a straight- 
forward solution of several interesting problems. Two geode- 
sic projections, the azimuthal equidistant projection and the 
Cassini-Soldner projection, are simple to write and their do- 
main of applicability is not artificially restricted, as would be 
the case, for example, if the se ries expansion for the Cassini- 
Soldner projection were used (Snvderl ll987l §13); the scales 
for these projections are simply given in terms of m\i and 
M12. Several other problems can be readily tackled with this 
library, e.g., solving other ellipsoidal trigonometry problems 
and finding the median line and other maritime boundaries. 
These and other problems are explored in lKarnevI d201 ll) . The 
web page http://geographiclib.sf.net/geod.html provides addi- 
tional information, including the Maxima (2009) code used to 
carry out the Taylor expansions and a JavaScript implemen- 
tation which allows geodesic problems to be solved on many 
portable devices. 



Acknowledgments 

I would like to thank Rod Deakin, John Nolton, Peter Os- 
borne, and the referees of this paper for their helpful com- 
ments. 



12 



References 

E. Beltrami, 1865, Risoluzione del problema: Riportare i punti di 
una superfi.de sopra un piano in modo die le linee geodetiche 
vengano rappresentate da linee rette, Annali Mat. Pura App., 7, 
185-204, |http://books.googlexom/boo ks?id=dfg EAAAAYAAJ&] 
|pg=PA185| 

F. W. Bessel, 1825, fiber die Berechnung der geographis- 
chen Ldngen und Breiten aus geoddtischen Vermessungen, 
Astron. Nachr., 4(86), 241-254, http://adsabs.harvard.edu/abs/ 

1825 AN 4..241B1 translated into English by C. F. F. Karney 

and R. E. Deakin as The calculation of longitude and latitude from 
geodesic measurements, Astron. Nachr. 331(8), 852-861 (2010), 
Idoi: 10. 1002/asna.20101 13521 E-print jarXiv:0908. 18241 

B. R. Bowring, 1997, The central projection of the spheroid and sur- 
face lines, Survey Review, 34(265), 163-173. 

L. M. Bugayevskiy and J. P. Snyder, 1995, Map Projections: A Ref- 
erence Manual (Taylor & Francis, London), http://www.worldcat. 
|org/oclc/31737484| 

E. B. Christoffel, 1868, Allgemeine Theorie der geoddtischen 
Dreiecke, Math. Abhand. Konig. Akad. der Wiss. zu Berlin, 8, 
119-176, |http://books.google.com/books?id=EEtFAAAAcAAJ&| 
|pg=PA119| 

A. C. Clairaut, 1735, Determination geometrique de la perpendicu- 
laire a la meridienne tracee par M. Cassini, Mem. de l'Acad. Roy. 
des Sciences de Paris 1733, pp. 406-416, http://books. google. 
|com/books?id=GOAEAA AAQAAJ&pg=PA406P 

C. W. Clenshaw, 1955, A note on the summation of Chebyshev series, 
Math. Tables Aids Comput., 9(51), 1 1 8-1 20, [http://www,jstor,org/| 
stable/2002068. 

J. S. Danielsen, 1989, The area under the geodesic, Survey Review, 
30(232), 61-66. 

C. F. Gauss, 1828, Disquisitiones Generates circa Superficies 
Curvas (Dieterich, Gbttingen), http://books.google.com/books? 
id=bXOAAAAAMAAJ, translated into English by J. C. Morehead 
and A. M. Hiltebeitel as General Investigations of Curved Sur- 
faces of 1827 and 1825 (Princeton Univ. Lib., 1902),| http://books. 
google.com/books ?id=a 1 wT JR3kHwUC . 

F. R. Helmert, 1880, Die Mathematischen und Physikalischen Theo- 
rieen der Hdheren Geoddsie, volume 1 (Teubner, Leipzig), http:// 
books. google. com/books ?id=qt2CAAAAIAAJ, translated into 
English by Aeronautical Chart and Information Center (St. Louis, 
1964) as Mathematical and Physical Theories of Higher Geodesy, 
Part 1, http://geographiclib.sf.net/geodesic-papers/helmert80-en. 
html. 

C. G. J. Jacobi, 1891, Uber die Curve, welche alle von einem 
Punkte ausgehenden geoddtischen Linien eines Rotationsellip- 
soides beriihrt, in K. T. W. Weierstrass, editor, Gesammelte 
Werke, volume 7, pp. 72-87 (Reimer, Berlin), op. post., com- 



pleted by F. H. A. Wangerin, http://books. google. com/books ?id=_ 
|09tAAAAMAAJ&pg=PA72] 

C. F. F. Karney, 201 1, Geodesies on an ellipsoid of revolution, Tech- 
nical report, SRI International, E-print arXiv: 1 10 2.1215vT] 

— , 2012, GeographicLib, version 1.20, http://geographiclib.sf.net. 

A. M. Legendre, 1806, Analyse des triangles traces sur la surface 
d'un spheroide, Mem. de lTnst. Nat. de France, 1st sem., pp. 130- 
161, [http://books. google.com/books?id=-dOEAAAAQAAJ&| 
|pg=PA130-lA4] 

I. G. Letoval'tsev, 1963, Generalization of the gnomonic projection 
for a spheroid and the principal geodetic problems involved in 
the alignment of surface routes, Geodesy and Aerophotography, 
5, 271-274, translation of Geodeziya i Aerofotos'emka 5, 61-68 
(1963). 

Maxima, 2009, A computer algebra system, version 5.20.1, |http://| 
maxima.sf.net. 

F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, ed- 
itors, 2010, NIST Handbook of Mathematical Functions (Cam- 
bridge Univ. Press), http://dlmf.nist.gov. 

B. Oriani, 1806, Elementi di trigonometria sferoidica, Pt. 1, 
Mem. dellTst. Naz. Ital., 1(1), 118-198, http://books.google.com/ 
|books?id=SydFAAAAc AAJ&pg=PA 1 1 8] 

— , 1808, Elementi di trigonometria sferoidica, Pt. 2, Mem. 
dellTst. Naz. Ital., 2(1), 1-58, http://www.archive.org/stream/ 
memoriedellistit2 1 isti#page/ 1 . 

— , 1810, Elementi di trigonometria sferoidica, Pt. 3, Mem. 
dellTst. Naz. Ital., 2(2), 1-58, http://www.archive.org/stream/] 
memoriedellistit22isti#page/l . 

R. H. Rapp, 1993, Geometric geodesy, part II, Technical report, Ohio 
State Univ., http://hdl.handle.net/181 1/24409. 

J. P. Snyder, 1987, Map projection — a working manual, Profes- 
sional Paper 1395, U.S. Geological Survey, http://pubs.er.usgs. 
gov/publication/pp 1395 

H. Vermeille, 2002, Direct transformation from geocentric coordi- 
nates to geodetic coordinates, J. Geod., 76(9), 451-454, doi:10. 
1 1007/s00190- 002-0273^6] 

T. Vincenty, 1975a, Direct and inverse solutions of geodesies on the 
ellipsoid with application of nested equations, Survey Review, 
23(176), 88-93, addendum: Survey Review 23(180), 294 (1976), 
http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf. 

— , 1975b, Geodetic inverse solution between antipodal points, 
unpublished report dated Aug. 28, http://geographiclib.sf.net/ 
geodesic-papers/vincenty75b.pdf. 

P. Wessel and W. H. F. Smith, 2010, Generic mapping tools, 4.5.5, 
http://gmt.soest.hawaii.edu/. 

R. Williams, 1997, Gnomonic projection of the surface of an ellip- 
soid, J. Nav., 50(2), 314-320, d oi:10.1017/S037 3463300023936 . 



